package org.systemsbiology.rface.hadoop.math.em;

import org.apache.commons.math.util.FastMath;
import org.systemsbiology.util.KahanSummation;
import org.systemsbiology.util.Summator;

/**
 * Created by IntelliJ IDEA.
 * User: anorberg
 * Date: 10/6/11
 * Time: 2:45 PM
 *
 * A representation of a Gaussian kernel density estimator as a mean in an E-M algorithm.
 */
public class GaussMean1D implements EstimatedMean1D{
    private double mean;
    private double variance;

    public GaussMean1D(double m0, double v0){
        mean = m0;
        variance = v0;
    }

    public double rawWeight(double value){
        double diff = value - mean;
        double dsquare = diff * diff;
        double gauss = FastMath.exp(-(dsquare / (2 * variance)));
        double norm = 1.0 / FastMath.sqrt(2 * Math.PI * variance);
        return gauss * norm;
    }

    public double revise(Iterable<? extends PointWeightViewer1D> points){
        Summator total = new KahanSummation();
        Summator weight = new KahanSummation();
        for(PointWeightViewer1D point: points){
            total.add(point.getPosition() * point.getWeight());
            weight.add(point.getWeight());
        }
        double net = total.total();
        double totalWeight = weight.total();
        double newMean = net/totalWeight;
        Summator varianceSum = new KahanSummation();
        for(PointWeightViewer1D point: points){
            double diff = newMean - point.getPosition();
            varianceSum.add(point.getWeight() * diff * diff);
        }
        double newVariance = varianceSum.total() / totalWeight;
        double ret = FastMath.abs(mean - newMean);
        ret += FastMath.abs(variance - newVariance);
        mean = newMean;
        variance = newVariance;
        return ret;
    }

    public double getMean(){
        return mean;
    }

    public double getVariance(){
        return variance;
    }
}
